{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solver comparison\n",
    "\n",
    "In this tutorial, we demonstrate how to compare geostatistical solvers using GeoStats.jl. This feature addresses a recurrent, sometimes underestimated, issue in the geostatistics community: the absence of scientific methods and software for rigorous selection of solvers and parameter sets.\n",
    "\n",
    "As we will see, the ability to manipulate solvers as first-class objects can save us from exhausting, manual trial and error, and can definitely help us make better use of our time on a critical project with many time constraints.\n",
    "\n",
    "Before we proceed, please install the following packages:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mPackage GeoStats is already installed\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mPackage InverseDistanceWeighting is already installed\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mPackage Plots is already installed\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mPackage PyPlot is already installed\n",
      "\u001b[39m"
     ]
    }
   ],
   "source": [
    "for pkg in [\"GeoStats\", \"InverseDistanceWeighting\", \"Plots\", \"PyPlot\"]\n",
    "    Pkg.add(pkg)\n",
    "end\n",
    "\n",
    "# make sure this tutorial is reproducible\n",
    "srand(2017);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem definition\n",
    "\n",
    "We start by creating a simple 2D estimation problem based on a given data set:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2D EstimationProblem\n",
       "  data:      500×3 GeoDataFrame (x and y)\n",
       "  domain:    100×100 RegularGrid{Float64,2}\n",
       "  variables: permeability (Float64)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using GeoStats\n",
    "\n",
    "geodata = readtable(\"data/permeability.csv\", coordnames=[:x,:y])\n",
    "domain = bounding_grid(geodata, [100,100])\n",
    "problem = EstimationProblem(geodata, domain, :permeability)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's suppose that we want to compare Kriging and inverse distance weighting (a.k.a. IDW) on this problem. Suppose that we are also interested in comparing different variogram models in Kriging. We define the solvers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{GeoStatsBase.AbstractEstimationSolver,1}:\n",
       " GeoStats.Kriging                      \n",
       " GeoStats.Kriging                      \n",
       " InverseDistanceWeighting.InvDistWeight"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using InverseDistanceWeighting\n",
    "\n",
    "solver₁ = Kriging(:permeability => @NT(variogram=ExponentialVariogram(range=40.)))\n",
    "solver₂ = Kriging(:permeability => @NT(variogram=SphericalVariogram(range=40.)))\n",
    "solver₃ = InvDistWeight()\n",
    "\n",
    "solvers = [solver₁, solver₂, solver₃]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and proceed to the comparisons.\n",
    "\n",
    "## Visual comparison\n",
    "\n",
    "Visual comparisons can be useful for quickly pruning solver/parameter combinations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Plots; pyplot(size=(1000,1000))\n",
    "\n",
    "compare(solvers, problem, VisualComparison())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cross-validation\n",
    "\n",
    "For quantitative comparisons, k-fold cross validation is a good option. We can plot the distribution of errors for each estimation solver and decide which is the best for the problem at hand based on summary statistics like the absolute mean validation error:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results = compare(solvers, problem, CrossValidation(10))\n",
    "\n",
    "# get validation errors for permeability\n",
    "errors = results[:permeability]\n",
    "\n",
    "# plot error distribution for each solver\n",
    "plt = plot(size=(1000,400), layout=(3,1), link=:x)\n",
    "\n",
    "histogram!(plt[1], errors[1], bins=50, color=:blue, alpha=.5, label=\"Kriging + exponential variogram\")\n",
    "vline!(plt[1], [mean(errors[1])], color=:black, alpha=1., lw=2, ls=:dash, label=\"mean\")\n",
    "\n",
    "histogram!(plt[2], errors[2], bins=50, color=:green, alpha=.5, label=\"Kriging + spherical variogram\")\n",
    "vline!(plt[2], [mean(errors[2])], color=:black, alpha=1., lw=2, ls=:dash, label=\"mean\")\n",
    "\n",
    "histogram!(plt[3], errors[3], bins=50, color=:red, alpha=.5, label=\"Inverse distance weighting\", xlab=\"validation error\")\n",
    "vline!(plt[3], [mean(errors[3])], color=:black, alpha=1., lw=2, ls=:dash, label=\"mean\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We verify the unbiasedness property of Kriging estimators (mean of errors is zero). We also stress that choosing a solver simply based on the mean absolute validation error is not a good idea. Although in this example, all the solvers present mean of errors around zero, inverse distance weighting is overestimating permeability values at many locations of the domain (see the mode in the right side of the distribution)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.1-pre",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
